#include <iostream>
#include <assert.h>
#include <Defines.h>
#include <TimeDiff.h>
#include <QTSKernel.h>
#include <MallocAlign.h>
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
using namespace std;
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
const  U32 bd64rinc[] = { 0, 0, 2, 2, 0, 0, 2, 2, 1, 1, 3, 3, 1, 1, 3, 3 };
const  U32 bd64cinc[] = { 3, 1, 3, 1, 2, 0, 2, 0, 3, 1, 3, 1, 2, 0, 2, 0 };
const  U32 ad64rinc[] = { 3, 1, 3, 1, 2, 0, 2, 0, 3, 1, 3, 1, 2, 0, 2, 0 };
const  U32 ad64cinc[] = { 0, 0, 2, 2, 0, 0, 2, 2, 1, 1, 3, 3, 1, 1, 3, 3 };
const  U32 bd04rinc[] = { 0,  0,  1,  1 };
const  U32 bd04cinc[] = { 1,  0,  1,  0 };
const  U32 ad04rinc[] = { 1,  0,  1,  0 };
const  U32 ad04cinc[] = { 0,  0,  1,  1 };
const  U32 od04rinc[] = { 0,  1,  1,  0 };
const  U32 od04cinc[] = { 0,  1,  0,  1 };
const char od04vrgn[] = { OD, OD, BD, AD };
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
class QTSArgs {
public:
   
   U32 argIdx;
   U32 m;
   U32 row;
   U32 col;
   U32 hgt;
   U32 rootIdx;
   char region;
   bool empty;
   
   const U64 * _nodes;

   const double * _x;
   double * _y;
};
class QTSThreadArgs {
public:
   
   U32 tidx;
   U32 numArgs;
   vector< QTSArgs > args;
};
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#ifdef __QTS__
U32 Pos2Row( const U64 _pos, const U32 _m, const U32 _n ) {
   U32 r = 0;
   U32 c = 0;
   Pos2RowColQTS( r, c, _m, _n, _pos );
   return r + 1;
}
U32 Pos2Col( const U64 _pos, const U32 _m, const U32 _n ) {
   U32 r = 0;
   U32 c = 0;
   Pos2RowColQTS( r, c, _m, _n, _pos );
   return c + 1;
}
#else
U32 Pos2Row( const U64 _pos, const U32 _m, const U32 _n ) {
   U64 pos = _pos;
   U32 row = 0;
   const U32 n = _n - 1;
   while( pos > n ) {
      pos -= _n;
      row++;
   }
   return 1 + row;
}
U32 Pos2Col( const U64 _pos, const U32 _m, const U32 _n ) {
   return 1 + _pos % _n;
}
#endif
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
struct ThreadRegions {
   vector< U32 > rmin;
   vector< U32 > rmax;
   vector< U32 > cmin;
   vector< U32 > cmax;
};
void AssignWorkToQTSThreads( const U32 numThreads, vector< QTSThreadArgs > & targs, vector< QTSArgs > & vargs, const U32 numRegions, const U32 m2 ) {
   
   U32 argIdx = 0;
   for( U32 rbi = 0; rbi < numRegions; rbi++ ) {
      for( U32 cbi = 0; cbi < numRegions; cbi++ ) {
         
         U32 tidx = rbi % numThreads;
         
         targs[ tidx ].args.push_back( vargs[ argIdx ] );
         targs[ tidx ].tidx = tidx;
         targs[ tidx ].numArgs++;
         
         argIdx++;
      }
   }
}
U32 FindQTThreadRegions( ThreadRegions & rgns, const S32 numRegions, const U32 dim ) {

   const U32 inc = dim / numRegions;
   const U32 tot = numRegions*numRegions;
   
   vector< bool > assigned( tot, false );

   U32 numAssigned = 0;

   rgns.rmin.clear();
   rgns.rmax.clear();
   rgns.cmin.clear();
   rgns.cmax.clear();
   
   for( U32 idx = 0; idx < tot; idx++ ) {
      rgns.rmin.push_back( 0 );
      rgns.rmax.push_back( 0 );
      rgns.cmin.push_back( 0 );
      rgns.cmax.push_back( 0 );
   }
   
   for( S32 rbi = 0; rbi < numRegions; rbi++ ) {
      for( S32 cbi = 0; cbi < numRegions; cbi++ ) {
         U32 rmin = ( rbi + 0 ) * inc - 0;
         U32 rmax = ( rbi + 1 ) * inc - 1;
         U32 cmin = ( cbi + 0 ) * inc - 0;
         U32 cmax = ( cbi + 1 ) * inc - 1;
         
         rgns.rmin[ numAssigned ] = rmin;
         rgns.rmax[ numAssigned ] = rmax;
         rgns.cmin[ numAssigned ] = cmin;
         rgns.cmax[ numAssigned ] = cmax;
         
         numAssigned++;
      }
   }
   
   assert( numAssigned == tot );
   return tot;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void ExploreQTSDagForRoots( QTSArgs * pargs, vector< QTSArgs > & vargs, const ThreadRegions & threadRegions, const U32 totNumRegions ) {
   
   DEBUG_PRINT( "ExploreQTSDagForRoots()" );
   const U32 numLevelsDown = LogXCeil( totNumRegions, 4 );
   const U32 tgtLevel = pargs->hgt - numLevelsDown;
   
   DEBUG_PRINT( "... hgt = " << pargs->hgt );
   DEBUG_PRINT( "... numLevelsDown = " << numLevelsDown );
   DEBUG_PRINT( "... tgtLevel = " << tgtLevel );
   
   if( numLevelsDown == 0 ) {
      // -- NOTE -- 
      // return the root node; but the default argument is already be populated for the root node:
      // thus immediate return
      return;
   } 
   
   vector< bool > found( totNumRegions, false );
   
   const U32 _hgt = pargs->hgt;
   const U32 _m   = pargs->m;
   const U32 rootIdx = pargs->rootIdx;
   const U32 * nodes = ( const U32 * ) pargs->_nodes;
      
   const double * _x = pargs->_x;
   double * _y = pargs->_y;
   
   const U32 numAlloc = 40;

   U32    ptrStack[ numAlloc ];
   U32  adptrStack[ numAlloc ];
   U32  bdptrStack[ numAlloc ];
   
   char   lvlStack[ numAlloc ];
   char adlvlStack[ numAlloc ];
   char bdlvlStack[ numAlloc ];
   
   U32    rowStack[ numAlloc ];
   U32    colStack[ numAlloc ];
   U32    adrStack[ numAlloc ];
   U32    adcStack[ numAlloc ];
   U32    bdrStack[ numAlloc ];
   U32    bdcStack[ numAlloc ];
   
   const U32 pmin = 0x80000000;
   const U32 pmsk = pmin - 1;

   U32 lvl;
   U32 path;
   U32 r;
   U32 c;
   
   U32 sp;   // stack-ptr
   U32 adsp; // above-diag stack-ptr
   U32 bdsp; // below-diag stack-ptr

   if( pargs->region & OD ) {
      sp   = 1;
      adsp = 0;
      bdsp = 0;
      lvlStack   [ 0 ]   = _hgt;
      ptrStack   [ 0 ]   = rootIdx;
      rowStack   [ 0 ]   = pargs->row;
      colStack   [ 0 ]   = pargs->col;
   }
   else {
      assert( false );
   }
   
   while( sp ) {
      
      while( sp >= bdsp && sp >= adsp && sp ) {

         sp--;
         const LSU2X * p = ( const LSU2X * ) &nodes[ ptrStack[ sp ] ];

         lvl = lvlStack[ sp ];
         r = rowStack[ sp ];
         c = colStack[ sp ];
         
         if( lvl <= tgtLevel ) {
            // find tidx
            U32 tidx = 0;
            while( tidx < totNumRegions ) {
               if( r >= threadRegions.rmin[ tidx ] && r <= threadRegions.rmax[ tidx ] ) {
                  if( c >= threadRegions.cmin[ tidx ] && c <= threadRegions.cmax[ tidx ] ) {
                     break;
                  }
               }
               tidx++;
            }
            assert( tidx < totNumRegions );
            
            // store this node as a regional root node
            vargs[ tidx ].argIdx = tidx;
            vargs[ tidx ].m    = _m >> ( _hgt - lvl );
            vargs[ tidx ].row  = r;
            vargs[ tidx ].col  = c;
            vargs[ tidx ].hgt  = lvl;
            vargs[ tidx ].rootIdx = ptrStack[ sp ];
            vargs[ tidx ].region  = OD;
            vargs[ tidx ]._nodes  = pargs->_nodes;
            vargs[ tidx ]._x = _x;
            vargs[ tidx ]._y = _y;
            vargs[ tidx ].empty = false;
            found[ tidx ] = true;
            DEBUG_PRINT( "... (OD) found thread region " << SETDEC( 3 ) << tidx << ": (" << SETDEC( 2 ) << r << ", " << SETDEC( 2 ) << c << ") "<< SETHEX( 16 ) << p->l.l << " " << SETHEX( 16 ) << p->r.l );
         }
         else if( lvl ) {
               
            if( p->l.s.lo & pmin ) {

               // path compaction...
               char rcrgn = OD;

               register U32 rx = r;
               register U32 cx = c;

               path  = p->l.s.hi;
               
               while( path > 2 ) {
                  lvl--;
                  U32 rpath = path & 3;
                  
                  if( rcrgn & OD ) {
                     const  U32 rinc[] = {  0,  1,  1,  0 };
                     const  U32 cinc[] = {  0,  1,  0,  1 };
                     const char vrgn[] = { OD, OD, BD, AD };

                     #ifdef LOCAL_DEBUG_ENABLED
                     U32 ppr = rinc[ rpath ] << ( lvl );
                     U32 ppc = cinc[ rpath ] << ( lvl );
                     #endif
                     DEBUG_PRINT( "... ... rpath = " << rpath << ", rinc = " << ppr << ", cinc = " << ppc << ", r = " << r << ", c = " << c << ", rx = " << rx << ", cx = " << cx );

                     rx += rinc[ rpath ] << ( lvl );
                     cx += cinc[ rpath ] << ( lvl );
                     rcrgn = vrgn[ rpath ];
                  }
                  else if( rcrgn & AD ) {
                     const  U32 rinc[] = {  1,  0,  1,  0 };
                     const  U32 cinc[] = {  0,  0,  1,  1 };
                     
                     #ifdef LOCAL_DEBUG_ENABLED
                     U32 ppr = rinc[ rpath ] << ( lvl );
                     U32 ppc = cinc[ rpath ] << ( lvl );
                     #endif
                     DEBUG_PRINT( "... ... rpath = " << rpath << ", rinc = " << ppr << ", cinc = " << ppc << ", r = " << r << ", c = " << c << ", rx = " << rx << ", cx = " << cx );

                     rx += rinc[ rpath ] << ( lvl );
                     cx += cinc[ rpath ] << ( lvl );
                  }
                  else {
                     const  U32 rinc[] = {  0,  0,  1,  1 };
                     const  U32 cinc[] = {  1,  0,  1,  0 };

                     #ifdef LOCAL_DEBUG_ENABLED
                     U32 ppr = rinc[ rpath ] << ( lvl );
                     U32 ppc = cinc[ rpath ] << ( lvl );
                     #endif
                     DEBUG_PRINT( "... ... rpath = " << rpath << ", rinc = " << ppr << ", cinc = " << ppc << ", r = " << r << ", c = " << c << ", rx = " << rx << ", cx = " << cx );

                     rx += rinc[ rpath ] << ( lvl );
                     cx += cinc[ rpath ] << ( lvl );
                  }
                  path >>= 2;
               }
               
               if( rcrgn & OD ) {
                  rowStack [ sp ] = rx;
                  colStack [ sp ] = cx;
                  ptrStack [ sp ] = p->l.s.lo & pmsk;
                  lvlStack [ sp ] = lvl;
                  sp++;
               }
               else if( rcrgn & AD ) {
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->l.s.lo & pmsk;
                  adlvlStack [ adsp ] = lvl;
                  adsp++;
               }
               else if( rcrgn & BD ) {
                  bdrStack  [ bdsp ] = rx;
                  bdcStack  [ bdsp ] = cx;
                  bdptrStack[ bdsp ] = p->l.s.lo & pmsk;
                  bdlvlStack[ bdsp ] = lvl;
                  bdsp++;
               }
            }
            else {
               
               U32 rclvl = lvl - 1;
               U32 sa = _hgt - rclvl;               
               const U32 inc = _m >> sa;
               
               if( p->l.s.lo ) {
                  // on diag: A_11
                  U32 rx = r;
                  U32 cx = c;
                  rowStack [ sp ] = rx;
                  colStack [ sp ] = cx;
                  ptrStack [ sp ] = p->l.s.lo;
                  lvlStack [ sp ] = rclvl;
                  sp++;
               }
               if( p->l.s.hi ) {
                  // on diag: A_12
                  U32 rx = r + inc;
                  U32 cx = c + inc;
                  rowStack [ sp ] = rx;
                  colStack [ sp ] = cx;
                  ptrStack [ sp ] = p->l.s.hi;
                  lvlStack [ sp ] = rclvl;
                  sp++;                  
               }
               if( p->r.s.lo ) {
                  // below diag: A_21
                  U32 rx = r + inc;
                  U32 cx = c;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->r.s.lo;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
               if( p->r.s.hi ) {
                  // above diag: A_12
                  U32 rx = r;
                  U32 cx = c + inc;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->r.s.hi;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
            }
         }
         else {
            // leaf
            assert( false );
         }
      }
      while( bdsp ) {

         bdsp--;
         const LSU2X * p = ( const LSU2X * ) &nodes[ bdptrStack[ bdsp ] ];

         lvl = bdlvlStack[ bdsp ];
         r = bdrStack[ bdsp ];
         c = bdcStack[ bdsp ];

         if( lvl <= tgtLevel ) {
            // find tidx!!
            U32 tidx = 0;
            while( tidx < totNumRegions ) {
               if( r >= threadRegions.rmin[ tidx ] && r <= threadRegions.rmax[ tidx ] ) {
                  if( c >= threadRegions.cmin[ tidx ] && c <= threadRegions.cmax[ tidx ] ) {
                     break;
                  }
               }
               tidx++;
            }
            assert( tidx < totNumRegions );
            
            vargs[ tidx ].argIdx = tidx;
            vargs[ tidx ].m    = _m >> ( _hgt - lvl );
            vargs[ tidx ].row  = r;
            vargs[ tidx ].col  = c;
            vargs[ tidx ].hgt  = lvl;
            vargs[ tidx ].rootIdx = bdptrStack[ bdsp ];
            vargs[ tidx ].region  = BD;
            vargs[ tidx ]._nodes  = pargs->_nodes;
            vargs[ tidx ]._x = _x;
            vargs[ tidx ]._y = _y;
            vargs[ tidx ].empty = false;
            found[ tidx ] = true;
            DEBUG_PRINT( "... (BD) found thread region " << SETDEC( 3 ) << tidx << ": (" << SETDEC( 2 ) << r << ", " << SETDEC( 2 ) << c << ") "<< SETHEX( 16 ) << p->l.l << " " << SETHEX( 16 ) << p->r.l );
         }
         else if( lvl ) {
                        
            if( p->l.s.lo & pmin ) {
               
               // path compaction...
               // DEBUG_PRINT( "... entering a path compaction, from below diagonal" );

               U32 rx = r;
               U32 cx = c;

               path  = p->l.s.hi;
               
               while( path > 2 ) {
                  lvl--;
                  U32 rpath = path & 3;
                  
                  const  U32 rinc[] = {  0,  0,  1,  1 };
                  const  U32 cinc[] = {  1,  0,  1,  0 };

                  #ifdef LOCAL_DEBUG_ENABLED
                  U32 ppr = rinc[ rpath ] << lvl;
                  U32 ppc = cinc[ rpath ] << lvl;
                  #endif
                  DEBUG_PRINT( "... ... rpath = " << rpath << ", rinc = " << ppr << ", cinc = " << ppc << ", r = " << r << ", c = " << c << ", rx = " << rx << ", cx = " << cx );

                  rx += rinc[ rpath ] << lvl;
                  cx += cinc[ rpath ] << lvl;

                  path >>= 2;
               }
               
               bdrStack  [ bdsp ] = rx;
               bdcStack  [ bdsp ] = cx;
               bdptrStack[ bdsp ] = p->l.s.lo & pmsk;
               bdlvlStack[ bdsp ] = lvl;
               bdsp++;
            }
            else {

               U32 rclvl = lvl - 1;
               U32 sa = _hgt - rclvl;
               const U32 inc = _m >> sa;
               
               if( p->l.s.lo ) {
                  U32 rx = r;
                  U32 cx = c + inc;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->l.s.lo;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
               if( p->l.s.hi ) {
                  // on diag: A_12
                  U32 rx = r;
                  U32 cx = c;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->l.s.hi;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
               if( p->r.s.lo ) {
                  // below diag: A_21
                  U32 rx = r + inc;
                  U32 cx = c + inc;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->r.s.lo;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
               if( p->r.s.hi ) {
                  // above diag: A_12
                  U32 rx = r + inc;
                  U32 cx = c;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->r.s.hi;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
            }
         }
         else {
            // leaf
            assert( false );
         }
      }
      while( adsp ) {
         
         adsp--;
         const LSU2X * p = ( const LSU2X * ) &nodes[ adptrStack[ adsp ] ];

         lvl = adlvlStack[ adsp ];
         r = adrStack[ adsp ];
         c = adcStack[ adsp ];

         if( lvl <= tgtLevel ) {
            // find tidx!!
            U32 tidx = 0;
            while( tidx < totNumRegions ) {
               if( r >= threadRegions.rmin[ tidx ] && r <= threadRegions.rmax[ tidx ] ) {
                  if( c >= threadRegions.cmin[ tidx ] && c <= threadRegions.cmax[ tidx ] ) {
                     break;
                  }
               }
               tidx++;
            }
            assert( tidx < totNumRegions );
            
            vargs[ tidx ].argIdx = tidx;
            vargs[ tidx ].m    = _m >> ( _hgt - lvl );
            vargs[ tidx ].row  = r;
            vargs[ tidx ].col  = c;
            vargs[ tidx ].hgt  = lvl;
            vargs[ tidx ].rootIdx = adptrStack[ adsp ];
            vargs[ tidx ].region  = AD;
            vargs[ tidx ]._nodes  = pargs->_nodes;
            vargs[ tidx ]._x = _x;
            vargs[ tidx ]._y = _y;
            vargs[ tidx ].empty = false;
            found[ tidx ] = true;
            DEBUG_PRINT( "... (AD) found thread region " << SETDEC( 3 ) << tidx << ": (" << SETDEC( 2 ) << r << ", " << SETDEC( 2 ) << c << ") "<< SETHEX( 16 ) << p->l.l << " " << SETHEX( 16 ) << p->r.l );
         }
         else if( lvl ) {
            
            if( p->l.s.lo & pmin ) {

               // path compaction...
               U32 rx = r;
               U32 cx = c;
               
               path  = p->l.s.hi;
               
               while( path > 2 ) {
                  lvl--;
                  U32 rpath = path & 3;
                  
                  const  U32 rinc[] = {  1,  0,  1,  0 };
                  const  U32 cinc[] = {  0,  0,  1,  1 };
                  
                  rx += rinc[ rpath ] << ( lvl );
                  cx += cinc[ rpath ] << ( lvl );
                  
                  path >>= 2;
               }

               adrStack   [ adsp ] = rx;
               adcStack   [ adsp ] = cx;
               adptrStack [ adsp ] = p->l.s.lo & pmsk;
               adlvlStack [ adsp ] = lvl;
               adsp++;
               
            }
            else {

               U32 rclvl = lvl - 1;
               U32 sa = _hgt - rclvl;
               const U32 inc = _m >> sa;
               
               if( p->l.s.lo ) {

                  U32 rx = r + inc;
                  U32 cx = c;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->l.s.lo;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
               if( p->l.s.hi ) {

                  U32 rx = r;
                  U32 cx = c;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->l.s.hi;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
               if( p->r.s.lo ) {

                  U32 rx = r + inc;
                  U32 cx = c + inc;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->r.s.lo;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
               if( p->r.s.hi ) {

                  U32 rx = r;
                  U32 cx = c + inc;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->r.s.hi;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
            }
         }
         else {
            // leaf
            assert( false );
         }
      }
   }
   for( U32 tidx = 0; tidx < totNumRegions; tidx++ ) {
      if( !found[ tidx ] ) {
         vargs[ tidx ].empty = true;
      }
   }
}
void SpMVFlatKernelQTS( QTSArgs * pargs ) {
      
   DEBUG_PRINT( "SpMVFlatKernelQTS()" );
   
   const U32 _hgt = pargs->hgt;
   const U32 _m   = pargs->m;
   const U32 rootIdx = pargs->rootIdx;
   const U32 * nodes = ( const U32 * ) pargs->_nodes;
   
   #ifdef __SPMV__
   const double * _x = pargs->_x;
   double * _y = pargs->_y;
   #endif
      
   const U32 numAlloc = 40;

   U32    ptrStack[ numAlloc ];
   U32  adptrStack[ numAlloc ];
   U32  bdptrStack[ numAlloc ];
   
   char   lvlStack[ numAlloc ];
   char adlvlStack[ numAlloc ];
   char bdlvlStack[ numAlloc ];
   
   U32    rowStack[ numAlloc ];
   U32    colStack[ numAlloc ];
   U32    adrStack[ numAlloc ];
   U32    adcStack[ numAlloc ];
   U32    bdrStack[ numAlloc ];
   U32    bdcStack[ numAlloc ];
   
   const U32 pmin = 0x80000000;
   const U32 pmsk = pmin - 1;

   U32 lvl;
   U32 path;
   U32 r;
   U32 c;
   
   U32 sp   = 0; // stack-ptr
   U32 adsp = 0; // above-diag stack-ptr
   U32 bdsp = 0; // below-diag stack-ptr

   if( pargs->region & OD ) {
      
      sp = 1;

      lvlStack   [ 0 ]   = _hgt;
      ptrStack   [ 0 ]   = rootIdx;
      rowStack   [ 0 ]   = pargs->row;
      colStack   [ 0 ]   = pargs->col;
      DEBUG_PRINT( "... (OD) r = " << pargs->row << ", c = " << pargs->col );
   }
   else if( pargs->region & BD ) {

      bdsp = 1;

      bdlvlStack [ 0 ]   = _hgt;
      bdptrStack [ 0 ]   = rootIdx;
      bdrStack   [ 0 ]   = pargs->row;
      bdcStack   [ 0 ]   = pargs->col;
      DEBUG_PRINT( "... (BD) r = " << pargs->row << ", c = " << pargs->col );
   }
   else {

      adsp = 1;

      adlvlStack [ 0 ]   = _hgt;
      adptrStack [ 0 ]   = rootIdx;
      adrStack   [ 0 ]   = pargs->row;
      adcStack   [ 0 ]   = pargs->col;
      DEBUG_PRINT( "... (AD) r = " << pargs->row << ", c = " << pargs->col );
   } 
   
   while( sp || bdsp || adsp ) {
      
      while( sp >= bdsp && sp >= adsp && sp ) {

         sp--;
         const LSU2X * p = ( const LSU2X * ) &nodes[ ptrStack[ sp ] ];

         lvl = lvlStack[ sp ];
         r = rowStack[ sp ];
         c = colStack[ sp ];
         
         if( lvl ) {
               
            if( p->l.s.lo & pmin ) {

               // path compaction...
               register char rcrgn = OD;

               register U32 rx = r;
               register U32 cx = c;

               path  = p->l.s.hi;
               
               while( path > 2 ) {
                  lvl--;
                  U32 rpath = path & 3;
                  
                  if( rcrgn & OD ) {

                     rx += od04rinc[ rpath ] << ( lvl );
                     cx += od04cinc[ rpath ] << ( lvl );
                     rcrgn = od04vrgn[ rpath ];
                  }
                  else if( rcrgn & AD ) {
                     
                     rx += ad04rinc[ rpath ] << ( lvl );
                     cx += ad04cinc[ rpath ] << ( lvl );
                  }
                  else {

                     rx += bd04rinc[ rpath ] << ( lvl );
                     cx += bd04cinc[ rpath ] << ( lvl );
                  }
                  path >>= 2;
               }
               
               if( rcrgn & OD ) {
                  rowStack [ sp ] = rx;
                  colStack [ sp ] = cx;
                  ptrStack [ sp ] = p->l.s.lo & pmsk;
                  lvlStack [ sp ] = lvl;
                  sp++;
               }
               else if( rcrgn & AD ) {
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->l.s.lo & pmsk;
                  adlvlStack [ adsp ] = lvl;
                  adsp++;
               }
               else if( rcrgn & BD ) {
                  bdrStack  [ bdsp ] = rx;
                  bdcStack  [ bdsp ] = cx;
                  bdptrStack[ bdsp ] = p->l.s.lo & pmsk;
                  bdlvlStack[ bdsp ] = lvl;
                  bdsp++;
               }
            }
            else {
               
               U32 rclvl = lvl - 1;
               U32 sa = _hgt - rclvl;               
               const U32 inc = _m >> sa;
               
               if( p->l.s.lo ) {
                  // on diag: A_11
                  register U32 rx = r;
                  register U32 cx = c;
                  rowStack [ sp ] = rx;
                  colStack [ sp ] = cx;
                  ptrStack [ sp ] = p->l.s.lo;
                  lvlStack [ sp ] = rclvl;
                  sp++;
               }
               if( p->l.s.hi ) {
                  // on diag: A_12
                  register U32 rx = r + inc;
                  register U32 cx = c + inc;
                  rowStack [ sp ] = rx;
                  colStack [ sp ] = cx;
                  ptrStack [ sp ] = p->l.s.hi;
                  lvlStack [ sp ] = rclvl;
                  sp++;                  
               }
               if( p->r.s.lo ) {
                  // below diag: A_21
                  register U32 rx = r + inc;
                  register U32 cx = c;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->r.s.lo;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
               if( p->r.s.hi ) {
                  // above diag: A_12
                  register U32 rx = r;
                  register U32 cx = c + inc;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->r.s.hi;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
            }
         }
         else {
            
            #ifdef __SPMV__
            _y[r] += p->l.d * _x[c];
            #else
            DEBUG_PRINT( "" );
            std::cout << SETDEC( 0 ) << 1+r << " " << SETDEC(0) << 1+c << " " << SETSCI( 20 ) << p->l.d << std::endl;
            //std::cout << SETDEC( 0 ) << 1+r << " " << SETDEC(0) << 1+c << " " << SETDEC( 0 ) << p->l.d << std::endl;
            #endif
         }
      }
      while( bdsp ) {

         bdsp--;
         const LSU2X * p = ( const LSU2X * ) &nodes[ bdptrStack[ bdsp ] ];

         lvl = bdlvlStack[ bdsp ];
         r = bdrStack[ bdsp ];
         c = bdcStack[ bdsp ];

         if( lvl ) {
                        
            if( p->l.s.lo & pmin ) {
               
               // path compaction...
               register U32 rx = r;
               register U32 cx = c;

               path  = p->l.s.hi;
        
               while( path > 2 ) {
                  
                  lvl--;
                  U32 rpath = path & 3;

                  rx += bd04rinc[ rpath ] << lvl;
                  cx += bd04cinc[ rpath ] << lvl;

                  path >>= 2;
               }
               
               bdrStack  [ bdsp ] = rx;
               bdcStack  [ bdsp ] = cx;
               bdptrStack[ bdsp ] = p->l.s.lo & pmsk;
               bdlvlStack[ bdsp ] = lvl;
               bdsp++;
            }
            else {

               U32 rclvl = lvl - 1;
               U32 sa = _hgt - rclvl;
               const U32 inc = _m >> sa;
               
               if( p->l.s.lo ) {

                  register U32 rx = r;
                  register U32 cx = c + inc;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->l.s.lo;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
               if( p->l.s.hi ) {

                  register U32 rx = r;
                  register U32 cx = c;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->l.s.hi;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
               if( p->r.s.lo ) {

                  register U32 rx = r + inc;
                  register U32 cx = c + inc;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->r.s.lo;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
               if( p->r.s.hi ) {

                  register U32 rx = r + inc;
                  register U32 cx = c;
                  bdrStack   [ bdsp ] = rx;
                  bdcStack   [ bdsp ] = cx;
                  bdptrStack [ bdsp ] = p->r.s.hi;
                  bdlvlStack [ bdsp ] = rclvl;
                  bdsp++;
               }
            }
         }
         else {

            #ifdef __SPMV__
            _y[r] += p->l.d * _x[c];
            #else
            DEBUG_PRINT( "" );
            std::cout << SETDEC( 0 ) << 1+r << " " << SETDEC(0) << 1+c << " " << SETSCI( 20 ) << p->l.d << std::endl;
            //std::cout << SETDEC( 0 ) << 1+r << " " << SETDEC(0) << 1+c << " " << SETDEC( 0 ) << p->l.d << std::endl;
            #endif
         }
      }
      while( adsp ) {
         
         adsp--;
         const LSU2X * p = ( const LSU2X * ) &nodes[ adptrStack[ adsp ] ];

         lvl = adlvlStack[ adsp ];
         r = adrStack[ adsp ];
         c = adcStack[ adsp ];

         if( lvl ) {
            
            if( p->l.s.lo & pmin ) {

               // path compaction...
               register U32 rx = r;
               register U32 cx = c;
               
               path  = p->l.s.hi;
               
               while( path > 2 ) {

                  lvl--;
                  U32 rpath = path & 3;
                                    
                  rx += ad04rinc[ rpath ] << ( lvl );
                  cx += ad04cinc[ rpath ] << ( lvl );
                  
                  path >>= 2;
               }
               
               adrStack   [ adsp ] = rx;
               adcStack   [ adsp ] = cx;
               adptrStack [ adsp ] = p->l.s.lo & pmsk;
               adlvlStack [ adsp ] = lvl;
               adsp++;
               
            }
            else {

               U32 rclvl = lvl - 1;
               U32 sa = _hgt - rclvl;
               const U32 inc = _m >> sa;
               
               if( p->l.s.lo ) {

                  register U32 rx = r + inc;
                  register U32 cx = c;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->l.s.lo;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
               if( p->l.s.hi ) {

                  register U32 rx = r;
                  register U32 cx = c;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->l.s.hi;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
               if( p->r.s.lo ) {

                  register U32 rx = r + inc;
                  register U32 cx = c + inc;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->r.s.lo;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
               if( p->r.s.hi ) {

                  register U32 rx = r;
                  register U32 cx = c + inc;
                  adrStack   [ adsp ] = rx;
                  adcStack   [ adsp ] = cx;
                  adptrStack [ adsp ] = p->r.s.hi;
                  adlvlStack [ adsp ] = rclvl;
                  adsp++;
               }
            }
         }
         else {
            
            #ifdef __SPMV__
            _y[r] += p->l.d * _x[c];
            #else
            DEBUG_PRINT( "" );
            std::cout << SETDEC( 0 ) << 1+r << " " << SETDEC(0) << 1+c << " " << SETSCI( 20 ) << p->l.d << std::endl;
            //std::cout << SETDEC( 0 ) << 1+r << " " << SETDEC(0) << 1+c << " " << SETDEC( 0 ) << p->l.d << std::endl;
            #endif
         }
      }
   }
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#ifdef __USE_PTHREADS__
void * SpMVQTSThread( void * vpargs ) {

   QTSThreadArgs * pargs = ( QTSThreadArgs* ) vpargs;
   SetThreadAffinity( pargs->tidx );
   
   for( U32 argIdx = 0; argIdx < pargs->numArgs; argIdx++ ) {

      if( !pargs->args[ argIdx ].empty ) {
         
         SpMVFlatKernelQTS( &(pargs->args[ argIdx ]) );
      }   
   }
   return NULL;
}
void SpMVQTSThreadLauncher( const U32 numThreads, const U32 regionSize, vector< QTSThreadArgs > & targs, pthread_t * threads ) {
   
   DEBUG_PRINT( "SpMVQTSThreadLauncher()" );
   DEBUG_PRINT( "... numThreads = " << numThreads );
   
   for( U32 tidx = 0; tidx < numThreads; tidx++ ) {
      pthread_create( &threads[ tidx ], NULL, SpMVQTSThread, &targs[ tidx ] );
   }
   for( U32 tidx = 0; tidx < numThreads; tidx++ ) {
      pthread_join( threads[ tidx ], NULL );
   }
}
#endif
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void SpMVQTS( const MatrixProperties & mtxMetaData, const U64 * nodes, const double * x, double * y, const U32 numRuns, const U32 numThreads, timespec * clks ) {
   
   DEBUG_PRINT( "SpMVQTS()" );
   
   const U32 deg = 4;
   const U32 hgt = mtxMetaData.dagHgts[ 0 ];
   const U64 rootIdx = deg * mtxMetaData.rootNodeIdxs[ 0 ];
   
   const U32 m2 = RoundUpToPow2( mtxMetaData.m );
   const U32 n2 = RoundUpToPow2( mtxMetaData.n );
   assert( m2 == n2 );

   DEBUG_PRINT( "... m2 = " << m2 << ", m = " << mtxMetaData.m );
   DEBUG_PRINT( "... n2 = " << m2 << ", n = " << mtxMetaData.n );
   DEBUG_PRINT( "... hgt = " << hgt );
   DEBUG_PRINT( "... deg = " << deg );
   DEBUG_PRINT( "... rootIdx = " << rootIdx );
   DEBUG_PRINT( "... mtxMetaData.rootNodeIdxs[ 0 ] = " << mtxMetaData.rootNodeIdxs[ 0 ] );
   #ifndef __SPMV__
   std::cout << "%%MatrixMarket matrix coordinate real general" << std::endl;
   std::cout << mtxMetaData.m << " " << mtxMetaData.n << " " << mtxMetaData.nnz << std::endl;
   #endif

   QTSArgs args;
   args.hgt      = hgt - 1;
   args.empty    = false;
   args.m        = m2;
   args.row      = 0;
   args.col      = 0;
   args.region   = OD;
   args._nodes   = nodes;
   args.rootIdx  = rootIdx;
   args._x       = x;
   args._y       = y;
   
   ////////////////////////////////////////////////////////////////////////////////
   ////////////////////////////////////////////////////////////////////////////////
   // find the thread safe root nodes
   #ifdef __USE_PTHREADS__
   const U32 regionSize = 128;
      
   ThreadRegions threadRegions;
   U32 totRegions = FindQTThreadRegions( threadRegions, regionSize, m2 );
   
   QTSThreadArgs targs;
   targs.tidx = 0;
   targs.numArgs = 0;
   
   vector< QTSThreadArgs > threadArgs( numThreads, targs );

   pthread_t * threads = (pthread_t*) calloc( numThreads, sizeof( pthread_t ) );

   vector< QTSArgs > vargs( totRegions, args );
   DEBUG_PRINT( "vargs.size() = " << vargs.size() );
   
   ExploreQTSDagForRoots( &args, vargs, threadRegions, totRegions );
   AssignWorkToQTSThreads( numThreads, threadArgs, vargs, regionSize, m2 );
   #endif
   ////////////////////////////////////////////////////////////////////////////////
   ////////////////////////////////////////////////////////////////////////////////
      
   #ifdef __USE_PTHREADS__
   SpMVQTSThreadLauncher( numThreads, numThreads, threadArgs, threads );
   #else
   SpMVFlatKernelQTS( &args );
   #endif

   #ifdef __CHECK_SPMV_RESULT__
   PrintPTRVectorToMMFile( x, mtxMetaData.n, false, "x_qts.mtx" );
   PrintPTRVectorToMMFile( y, mtxMetaData.m, false, "y_qts.mtx" );
   #endif
   
   #ifdef __SPMV__
   U32 j = 0;
   // begin next spmv iter:
   while( j < numRuns ) {
      clock_gettime( CLOCK_MONOTONIC, clks ); clks++;

      #ifdef __USE_PTHREADS__
      SpMVQTSThreadLauncher( numThreads, numThreads, threadArgs, threads );
      #else
      SpMVFlatKernelQTS( &args );
      #endif

      clock_gettime( CLOCK_MONOTONIC, clks ); clks++;
      j++;
   }
   #endif
   #ifdef __USE_PTHREADS__
   free( threads );
   #endif
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
